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Abstract: Axions constitute a well-motivated dark matter candidate, and if PQ sym¬ 
metry breaking occurred after inflation, it should be possible to make a clean prediction 
for the relation between the axion mass and the axion dark matter density. We show that 
axion (or other global) string networks in 3D have a network density that depends loga¬ 
rithmically on the string separation-to-core ratio. This logarithm would be about 10 times 
larger in axion cosmology than what we can achieve in numerical simulations. We simulate 
axion production in the early Universe, finding that, for the separation-to-core ratios we 
can achieve, the changing density of the network has little impact on the axion production 
efficiency. 
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1 Introduction 

The axion [1-4] is a proposed particle, the angular excitation of a new “Peccei-Quinn” (PQ) 
field (f that would solve the strong CP problem [5-7] and which is also a very interesting 
dark matter candidate [8-10], thereby solving two puzzles with one mechanism. In this 
light, the study of the axion as a dark matter candidate is, in our view, well motivated. 
The axion model has one undetermined parameter, the vacuum value of y?, fa] the axion 
mass TTia scales as The value of fa also determines the amount of axion dark matter 
that would be produced in the early Universe, which means that (for a well-motivated 
initial condition we will describe) it should be possible to predict the axion mass from the 
dark matter density. To do so, we need to understand the efficiency of axion production 
in cosmology. 
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As we will soon explain, the PQ field’s evolution in the early Universe is complicated 
by the appearance of structures - cosmic strings - that may play a role in determining 
the production efficiency of axions around the QCD scale T ~ 1 GeV. Here we will study 
these structures and their role in axion production via numerical simulations. Ours is not 
the first study on this problem [11-27], but we believe it will also not be the last; we will 
demonstrate some significant challenges for such studies, one of which we are not yet able 
to overcome. Axion cosmic string networks are in a family called “global string networks,” 
and Martins and Shellard have argued [27, 28] that such networks are sensitive to the 
sizes of the string cores, which cannot be properly simulated numerically. This implies 
that numerical simulations will view string networks with very different properties - in 
particular, a much lower string density - than the ones that would really occur for physical 
axion field parameter values. We will present strong numerical evidence supporting this 
claim. In particular, we will demonstrate that the cosmic string density increases as we 
increase the log of the ratio of the network age to the string core size. This ratio is very 
large in cosmology, implying that physically relevant string networks may be an order of 
magnitude denser than those in numerical axion-production simulations. We will also point 
out another potential pitfall to numerical simulations, which can be overcome but must be 
monitored. 

We use the remainder of the introduction to review the properties of axions relevant 
to cosmology and the overall picture of axion production in the early Universe. From the 
point of view of cosmology, the relevant features of the axion are that it is a complex scalar 
field if with a symmetry-breaking potential^, 

- ^axion = df,ip*d^ip + ^ [2ip*ip - flf (1.1) 

with fa the vacuum expectation value and A the self-coupling, which we assume is 0(1). 
The Lagrangian has a U(l) “PQ” symmetry ip —yje*®, spontaneously broken when p takes 
on a vacuum value somewhere on the “vacuum manifold” 2p*p = f’f. The symmetry is 
also explicitly broken by an anomalous coupling to QCD and instanton effects, giving rise 
to an extra contribution to the potential, 

- >Caxion,eff.QCD = X^) [1 “ COS (arg p)] . (1.2) 

Here x(r) is the temperature-dependent topological susceptibility of QCD, which for our 
purposes is an input from the theory of QCD. Its vacuum value^ is x(0) — (76 MeV)^. 

The minimum of the axion potential is at y/2 p = with (T = fa and 9a = 0. Radial 
fluctuations a = fa + s, which we will call saxions, have a mass = Xfa', and angular 
excitations 9a = a/fa, axions, have a mass m\{T) = x{T)/fa- If fa ^QCD, as required 
by existing constraints [31-33], then there is a large hierarchy between these masses, and 
only the axion should play a role cosmologically for T < /„. 

^Our metric convention is [—and we use standard complex-field notation p = {pr -t ipi)l\/2, 
explaining some strange-looking factors of 2. 

^At lowest order in chiral perturbation theory, yfT = 0) ~ Ffrn^ , with iT the pion decay 

constant, m,r the pion mass, and m„, md the up and down quark masses [29]. Recent lattice determinations 
[30] find iriulmd ^ 0.45. 
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x{T) is known to be a strong function of T at high temperature T ^ Aqcd; varying 
roughly as [ 34 ]. Therefore as time progresses and temperature falls cos- 

mologically, the axion goes rapidly from being effectively massless, mat <C 1 , to massive, 
mat S> 1. The dynamics around mat ~ tt are quite nontrivial; but once mat S> 1, axion 
fluctuations will be small and the axion number an adiabatic invariant, so the dynamics 
around mat ~ tt determine the efficiency of axion production and therefore the amount of 
axion dark matter. 

It seems likely that PQ symmetry is restored near the end of or after inflation. In 
particular, PQ symmetry is generically restored during high-scale inflation (if the Hubble 
parameter during inflation is large, H > fa) [35]^. Thermal symmetry restoration after 
inflation would have occurred if the Universe reheated to a temperature T > fa- If either 
occurred, then we know the initial conditions for the axion field: 6a would start out un¬ 
correlated at causally-disconnected points, meaning essentially random initial conditions 
should apply. With the initial conditions known statistically, the axion cosmology model 
has only one unknown parameter, /„. If we can determine x(T) and work out the axion 
dynamics near mat ~ vr, then it should be possible to compute the relation between the 
value of fa and the resulting dark matter density. If we assume that the axion consti¬ 
tutes the dark matter in the Universe, the known dark matter density [38] should give a 
unique prediction for fa, and with it the axion mass. The axion mass is experimentally 
measurable, making this scenario testable. And a narrow search window would also be 
very valuable for the design of the most sensitive type of experiment, resonant microwave 
cavity detectors [39-42]. 

To finish setting the stage, we describe qualitatively how the fields evolve near mat ~ vr. 
At early times, ~ 0 and the Lagrangian Eq. (1.1) has a 1/(1) symmetry. This symmetry 
is locally broken by the phase choice Random initial conditions give rise to a 

network of topologically stable cosmic strings [43], which then evolve and untangle, entering 
a scaling solution in which they maintain a string length per volume of order 

7-^string _ f, /, 

V f^ ’ ^ 

where 7, Tstring; and V are the typical Lorentz gamma-factor of a moving string, the 
physical length of string, and the physical volume under consideration, respectively, and 
t is the age of the Universe. ^ is an order-1 parameter describing the network evolution. 
After the axion mass becomes relevant, the potential has a single global minimum^, and 
there are no true topological structures. However, the strings remain metastable, and there 
are metastable domain walls associated with the field’s phase varying by 27r from one side 
of the wall to the other; exactly one wall ends on each string. The domain wall tension 
draws the strings together and accelerates the annihilation of the network [44]. 

®There may be ways around this argument [36, 37], but we believe that it is the generic expectation. 

"^One can also consider axions with multiple minima, so cosfarg <fi) in Eq. (1.2) becomes cos{N arg (p). 
However, it is difficult to avoid problems associated with stable domain walls in this model while simultane¬ 
ously solving the strong CP problem without fine tuning [15, 18, 19], so we will not consider this possibility 
here. 
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We will argue that the “scaling solution” for the strings is not quite as simple as 
Eq. (1.3) suggests; the string tension, the parameter and other properties of the string 
and string-wall networks should vary logarithmically with the ratio of the inter-string 
separation t/y/^ ~ t to the string core size ~ l/mg. In numerical held-theory simulations 
this ratio is bounded by the size of the lattice used to solve the field dynamics: rngt will 
not exceed about 10^ in 3D simulations or 10^ in 2D simulations. Physically, we want 
~ fa/H, with fa ~ 10^^ GeV and H ~ T'^/mpi ~ 10“^® GeV; so the physically 
relevant value for the ratio is more like rngt ~ 10^*^. For a fixed fa, the physical value 
of rUst gives strings with 10 times higher tension than in simulations, leading to a much 
denser string network which breaks up more slowly under the action of domain walls. The 
physical network evolution may be much more efficient at producing axions than recent 
simulations (including the ones we present), implying a smaller value of fa and heavier 
value of nia- 

The other danger is that even the meta-stability of the domain walls is dependent 
on the mass ratio malms- We will show below that, for m^/m^ > 1/39, the domain 
walls cease to be even metastable, and the topological structures abruptly collapse. This 
certainly should not happen cosmologically; but since ma rises rapidly with time, it is 
a challenge to make a numerical simulation with large enough that this condition is 
maintained until the network breaks up via the expected physics. Therefore, there is some 
danger that simulations will incorporate unphysical collapse of the string-wall network or 
must be stopped before the network has finished evolving. 

The next section reviews the physics of global cosmic strings and explains the rele¬ 
vance of ln{mst), using analytical arguments. Next we present numerical evidence that the 
string network grows denser with increasing In(mst). Then we study the network evolu¬ 
tion and axion production around mat ~ vr in more detail. Our study of the string-wall 
network’s collapse actually does not show evidence that the final axion density depends 
strongly on ln(m 5 t); but the dynamic range we can study is too narrow to make any brave 
extrapolations about what happens when the log is increased by another factor of 10. 

As an aside, we mention that besides the details of the axion dynamics around mat ~ vr, 
there are also uncertainties in the axion production efficiency from our incomplete knowl¬ 
edge of the temperature dependence of the topological susceptibility x{T)- The high- 
temperature behavior is only known at hrst order in perturbation theory [34], so the first 
unknown corrections are suppressed by 0{as)- Therefore the perturbative treatment may 
not be very reliable around T ~ 1 GeV where the interesting dynamics occurs. In this 
temperature range we only have model calculations [45, 46]; lattice calculations [47, 48] 
are currently available only at lower temperatures and/or in the quenched approximation. 
We will find that the axion production is not too sensitive to the exact value of x(T), but 
it would nevertheless be valuable to have a reliable lattice calculation for T in the range 
from 0.5 to 1.5 GeV. 
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2 String networks and string cores 


We begin with a lightning review of why the strings arising from Eq. (1.1) have logarithmi¬ 
cally large string tensions. We assume some familiarity with string defects; a reader who 
needs some background can look in [49]. We will make the classical field approximation 
throughout, which is an excellent approximation for the IR axion field dynamics since the 
mean occupancy is ~ fa/^^ ~ 10®®. 


2.1 Cosmic strings 

Consider a string lying along the z axis in polar coordinates. The field varies as (p = 
with cj) the azimuthal angle and /(r) a function obeying /(r) 0, 

f{r) 1 — 0(l/m^r^). The associated energy is 


^ — J = J dz J rdr J 


d(j) 




ffa 






(2.1) 

Far from the string’s core, / —)• 1 and all terms become negligible except for the (^-derivative 
term, which becomes /^/2r^. Therefore the energy density decays as 1/r^ as we move away 
from the string core. This is in contrast to “local” strings, where the U(l) symmetry is 
gauged (local) and a gauge field compensates for the (/)-derivative except in the string’s 
core. The local string case has received much more attention in the literature. 

The large-distance part of the string tension (energy per length) in Eq. (2.1) is ap¬ 
proximately 


T = 


E 

T 





( 2 . 2 ) 


where t is an IR length scale where this description breaks down, for instance the distance 
to the next string, £ ~ 1/H. This dependence on the distance to the next string indicates 
that there is a long-range force-per-length (gradient of string energy-per-length) between 
strings, with a strength of order d{E/L)/di = 

There are two important facts here. First, there are long-range interactions between 
strings, mediated by the (nearly) massless axion field. The force between strings of opposite 
winding sense is attractive, which helps them to hnd each other and annihilate. Though we 
have not shown it, the presence of a massless mode also helps accelerating or bent strings 
to radiate energy more efficiently than for the local string case. Both of these effects help 
the string network to annihilate more efficiently, leading to a much lower-density string 
network. 

Second, the attractive forces between strings, and the radiation of energy into long- 
wavelength axions, scale with in the same way as the string tension, but they are not 
enhanced by while the energy-per-length, or string tension, is. Therefore the ratio 

of tension to radiation/force effects is proportional to this log, but not to fa- To keep track 
of this difference, we will name this large logarithm k = ln(ms.^) ~ lii{fa/E[). 
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2.2 Scaling in 2D 

Let’s first see how k plays a role in string density for 2D networks, where sensitivity to this 
logarithm is generally accepted [11]. The 2D network is described by the 3D theory but 
enforcing that the field (/? does not vary along the z axis, ip = ip{x,y). In the x — y plane, 
the string becomes a monopole or point-charge, with ip{x, y) varying by ±27r as one goes 
around the charge, depending on whether the string has positive or negative orientation. 
The two orientations of strings act like two signs of charges. 

The analogy to electric charges turns out to be complete [50-54] . Outside of the string 
cores, we can write ip = and the equation of motion becomes 

= or d^d^^ea = o. (2.3) 

Except in string cores, we can define dual electric and magnetic variables 

^ , (2.4) 


obeying 

= fae^^^d^d^6a = 0 , = -2fad^d^9a = 0 , (2.5) 

which are the free-space Maxwell equations in 2-|-l dimensions. A surface in 2D is a loop, 
and the electric flux through the surface is 27r times the winding number around the loop. 



'C 


^ij Fidlj 



2vr/aUencl. 


showing that a string is the source for a flux of ±27r/a. 
The electrical attraction between two strings will be 


gig2 ^ _^ 27r/^ 
2'Kr r 


( 2 . 6 ) 


(2.7) 


as expected if each has an energy 7rf^ln[msr) as we found above^. The short-distance 
part of this energy, tt/^ In(msro) = vr/f k, should be interpreted as the mass of the charge, 
m = KT^fa- Varying k (the log of the ratio of separation to core length scales) is varying 
the charges’ mass at fixed charge magnitude. We want to argue that making the charges 
heavier will make them more non-relativistic, so they move and annihilate less efficiently 
and have a higher density. 

Suppose that at t = 0 there is a very dense starting ensemble of positive and nega¬ 
tive charges, which evolve under Hubble drag and their electromagnetic interactions. We 
can find how the density of charges will evolve by making parametric scaling estimates. 
The charges find each other under the influence of Coulomb attraction and annihilate off. 
Suppose that at time t the mean inter-charge separation is i, so the density of charges is 
and let us try to estimate i. On dimensional grounds we expect I ~ with n an 

exponent we want to determine. 

We see that the density of charges should fall, with 0(1) of the charges annihilating 
in an 0{t) amount of time. To do so, the charges have to move a distance ~ ^ in a time 

®The units may look strange, because an energy in 2D is an energy-per-length or tension in 3D units. 
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t, implying that v ~ ill. The typical force on a charge is T ~ 27rf^/£, so the kinetic 
energy a charge obtains is force times distance, e = F x d 2vr/2/^ X £ 27r/^. Equating 

with we find v ~ 2j^/K. Then £ ~ 2K~^^‘^t. The density of strings 

is n ~ £~‘^ ~ njAt'^. Therefore we estimate that the number density of strings should 
increase linearly with the log of the scale ratio k, and that the velocity should scale as 
1/y^, becoming nonrelativistic as the logarithm becomes large. 

We therefore have a robust argument that, in 2 space dimensions, the string density 
will scale with the log of the scale separation k = In(mst) as n ~ with small squared 
velocities (u^) ~ 1 /k. 


2.3 Scaling in 3D 


The argument in 2D clearly does not translate simply into 3D, since 3D strings move under 
tension as well as under mutual interactions. Nevertheless we expect strong, though not 
necessarily linear, ^-dependence in the string network density in 3D as well. The reason 
is that the long-range interactions of the strings provide rather efficient mechanisms for 
the strings to radiate energy and to annihilate against each other. This is in contrast 
to local (gauge) string networks, where there are no long-range interactions and only very 
inefficient radiation of gravitational waves. As a result, the density of global strings, ^ < 1.2 
at achievable k values, is an order of magnitude smaller than the value ^ ~ 13 found in 
numerical lattice field theory [27, 55] and Nambu-Goto [24, 56, 57] evolutions. 

However, while the interactions between strings and the radiation of energy from ac¬ 
celerating strings should scale as /f, the tension of strings should scale as f^K. In the 
large-K limit, the tension should dominate the mutual string interactions and radiation. 
Therefore, as k is increased global string networks should behave more like Nambu-Goto 
string networks. For very large values of k, the network density ^ should be a factor of 10 
larger than at currently achievable k values. 


Martins-Shellard modeled ((/c) 



«:=log(m^ Tg) 



/c=log(mgTo) 


Figure 1. Estimated k dependence of the string density ^ and velocity (v) according to the one-scale 
model of Martins and Shellard. 

Martins and Shellard attempted to model this by amending their 1-scale model for 
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string networks to include the radiation effects associated with global strings [28]. Fitting 
the single parameter describing radiation efficiency in their Eq.(4.2) to the value of ^ = 1.15 
at K = 6 from Fig. 3 below, and solving their equations describing the string density and 
velocity, we get the results plotted in Fig. 1. According to this model, the string density 
should vary nearly linearly with k for k ~ 5-10, and should grow by roughly a factor of 5 
in going from «; = 6to«; = 70. 

We will compare numerical simulations to these expectations in the next section, by 
considering string evolution without the “tilt” term, Eq. (1.2). 

3 Simulations of string networks 
3.1 Implementation 

To investigate these issues, we have implemented an MPI-parallelized code to evolve the 
equations of motion that follow from Eq. (1.1) and Eq. (1.2). We work in a radiation 
dominated FRW background, so the Hubble parameter® H = dtan/an is related to the 
time t as H = l/(2t). We work in comoving coordinates and introduce conformal time 

dr = —, so that r = 2i, (3.1) 

O'H V ®iro 

with to and a^o tke values of t and at some (arbitrary) hxed time. Similarly, we rescale 
masses by so that ruphys.dt = mconi.dT] in what follows, all masses are rriconi.- The 
temperature scales as T (x t~^, and the metric is g^u = In these 

units and writing space {i,j) and time (0) components and factors of r from the metric 
explicitly, the action, up to an irrelevant multiplicative constant, is 

d^x ^ — T^dr‘^*dr‘^ + T‘^diip*dnp (3.2) 

+ ^( 27 ?V-/a) +T^x('r)(l-cos(arg(/7))^ . 

We discretize this on a grid that is uniform in these coordinates, with spacing a and tem¬ 
poral spacing at = a/nt- We give a few more details about our numerical implementation 
in Appendix A; in particular we explain there how we determine the density of strings and 
walls, and what we believe is a new technique for establishing the string velocity directly 
from information in the field variables. 

Our numerical implementation is rather standard, except for two points. First, like 
many authors we replace 

X(T)[1-cos(arg(^)] ^ x(r)[l -/"Vrj , (3.3) 

with ifr = y/2 Re (/J. This form has the advantage of being analytic in held variables, so no 
trigonometric evaluations are required. It differs from the original form only in string cores, 
where we expect this term to have little impact; the symmetry-breaking term will always be 

®We write the Hubble parameter as an to avoid confusion with the lattice spacing a. 
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small, so it is only important because it applies over large regions of space where the field is 
near its minimum. We assume that x(t) grows as a power of r in the relevant time regime, 
so if x(^) oc then r^x(r) = /«two powers higher than the temperature 
dependence of x(T). We also scale out an overall factor of from the Lagrangian and 
rescale if ^ ip/fa- 

The other nonstandard change we make is to the potential term. Physically, we are 
interested in large values of = A/^. We cannot make the value larger than ~ 

1, since otherwise we encounter numerical artifacts (as discussed in Appendix A); but 
physically we want it to be as large as possible, since we expect the physical value of mg to 
be very large. Therefore we remove two powers of r from this term so that its r scaling is 
the same as the gradient terms, so m^a remains fixed in our simulations. In other words, 
we keep the string core size fixed in lattice units, rather than allowing it to grow smaller in 
comoving coordinates due to Hubble expansion. This change enlarges the dynamic range 
where our simulations can see cosmic string networks, making it easier for us to achieve 
“scaling.” 

3.2 Results for string-only networks 

As a first application, we investigate string networks without the symmetry-breaking term, 
that is, setting x(t) = 0 or = 0. In this case, no domain walls form, and the string 
network evolves until we terminate the simulation at r < L/ 2 . 

We have plotted the length (number) of strings, scaled by to account for system 
expansion, for 3D and 2D simulations in Fig. 2. The left plots are in terms of conformal 
time measured in units of the saxion mass, they show a rise, at first rapid and then 
more gradual, in the string density when normalized by a factor. To test the hypothesis 
that this rise represents logarithmic scaling in the separation-to-core ratio k ~ In(rms), 
we present a semilog plot on the right-hand side of the figure. Indeed, the nearly straight- 
line behavior is striking in both 2D and 3D. In each figure we have shown two or more 
choices of lattice spacing m^a, and we indicate the la statistical errors with upper and 
lower thin curves accompanying each best-value thick curve. The figure is also a check 
that our rather coarse lattice, anig = 1.5, is sufficient to see continuum behavior; we plot 
curves for anig = 1.5, for arug = 1.0, and (in 3D) for arug = 0.7, which fall on top of each 
other to within the small statistical errors. 

It is common in the cosmic string community to consider also the mean velocity, 
and Lorentz 7 -factor of the strings. It is more customary to normalize the network density 
in terms of the time rather than the conformal time; according to Eq. (3.1) this divides^ 
the string densities shown in Fig. 2 by 4. In addition, it is customary to track not the 
length of string network, but the energy in the network divided by the string tension. This 
is the same as integrating f 'jdl rather than f dl in determining the effective length of 
string in the network. In computing the mean velocity, and 7 -factor one also uses this 

^The factor is still rather than t, because one must also work in terms of physical rather than comoving 
lengths and volumes. 
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3D scaled string density 



3D scaled string density 





Figure 2. String length per volume, scaled by for 3D (upper) and 2D (lower) simulations, 
plotted against conformal time (left) and its log (right). Upper/lower curves are ±lcr statistical 
errors. 


normalization, so 


((^^) > {v^) > ( 7 )) 


J j dl X yv , v'^ 
J-fdl 



(3.4) 


We plot the normalized network density in 3D and 2D, using this normalization, in 
Fig. 3. The figures indicate that the 7 -weighted string density is also a logarithmic function 
of Tnig, and that the string density in 3D is over an order of magnitude smaller than the 
value ~ 13 for local string networks [56, 57]. 

We also plot the values of (u), (u^), and ( 7 ) in Fig. 4. The string velocity falls off at 
large rm<j in 2D as we expected, but in 3D it roughly approaches a constant, which is quite 
close to the value from Fig. 1. The figure shows that, while = 1.5 was a sufficient value 
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3D string ^-value 




Figure 3. String network density in 3D (left) and 2D (right), varying with the log of conformal 
time. Different curves refer to different lattice spacings rusa. 



2D Velocity, v^, 7 



Figure 4. String mean velocity, and 7 -factor in 3D (left) and 2D (right), varying with the log 

of conformal time. 


to obtain continuum-limit string densities (Figure 2 and 14), our algorithm for finding the 
string velocity (see Subsection A.2) is more sensitive to the size of the string core, and 
a high-quality determination of the string velocity and especially of its 7 -factor by this 
method requires m^a < 1 in 3D. We also see oscillations in string velocity at small nisT] 
we believe that these are string-core “breathing” modes induced by our initial conditions, 
which may fool our string-velocity method, since it assumes that the string core takes its 
unperturbed form. 

This section’s results show that the string density in 2D evolves exactly as we would 
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expect: the string density grows with k = In(rms), while the string {v^) falls as k~^. In 3D 
the string network density grows in a way which, for n values we have available, is roughly 
consistent with linear behavior. The string velocity in 3D is a very weak function of k and 
is near the value predicted by the one-scale model [28]. 


4 String-wall network evolution 
4.1 Scaling expectations 


Consider the axion held in conformal coordinates. As we discussed, the axion mass grows 
strongly as temperature drops, and therefore as conformal time progresses. We will desig¬ 
nate a special time tq, dehned as when 

ma(ro)ro = 1 (in terms of time, ma{to) = H{to )). (4.1) 


Roughly speaking, this is when the explicit symmetry breaking and axion mass start be¬ 
coming physically signihcant. We will parameterize the topological susceptibility near this 
time as x(T) oc T~'^. In our numerical work we will take n = 7. 

When ma{T)T ~ tt the dynamics are complex and must be solved numerically. But for 
n = 7, by T = 3 to we have ma{T)T ~ 420, which should provide enough held oscillations 
to force the string-wall network evolution to completion. By this time {dma/dT)/m‘^ <C 1, 
ensuring adiabatic evolution of the axion held. In our comoving and conformal coordinates, 
one then expects adiabatic evolution of the form e oc niafa/T'^ and® Uax — oc /^/r^. 

On dimensional grounds, in the radiation era and before the QCD transition when the 
number of radiation degrees of freedom changes, the axion number should be 

_ iVax _ Krofi 
V 

with K a constant. This constant determines the produced density of axions; since little 
entropy is produced between GeV temperatures and today, the axion-to-entropy ratio in 
the modern Universe is 



^ ^ nax(r = To) ^ KH{To)fl 
S s(To) ^ ^ 

where Tq is the temperature when t = tq and is the effective number of light degrees of 
freedom at T = Tq. Therefore, determining K is determining the key input for the axion 
abundance in the modern universe. 

Our goal is to determine this constant by evolving the axionic held, starting with a 
random space-varying phase at very early time, until the string network is gone and the 
axion huctuations are small, and then evaluating 


r d^k Ek 
Tq J (27r)3 Eh 


d^k 

(27r)® 



{oim + 


1 



(4.4) 


®The power may look strange. In the radiation epoch the conformal time is proportional to the 
expansion factor, so one would expect riax oc r“®. But the relation between time and conformal time means 
that a fixed particle number or energy appears to be growing as in conformal time. Also, we include 
a factor of the conformal-time rria in our scaling relation; and because of the time-scaling, a fixed rria in 
regular time is rua oc r in conformal time. Together these effects provide '7--3+2-1 _ 
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which should be independent of the final time r when we evaluate it, if that time comes 
after adiabatic behavior sets in. Eq. (4.4) can be evaluated by fast Fourier transform (EFT) 
methods, or by the method of Appendix B when the EFT is not available. 

As a baseline for the expected value of A', we can consider the angle-average of the 
misalignment mechanism; we evolve a spatially uniform initial condition for 6a, leading to 

d‘^9a = --drOa - sin( 6 »a) (4.5) 

r 

with rria = / Tq~^^ , and then take the average of the resulting axion number over values 

of the starting angle 6a- For our choice n = 7 we find K = 16.0255. 

4.2 String-wall network evolution in 2D 

We have evolved the lattice axion equations of motion for both 2D and 3D systems, using 
the coarsest lattice that gives continuum-like string behavior, mga = 1.5 (see App. A. 3). 
We consider a number of values for the dimensionless ratio rngTo. Depending on one’s 
perspective, a large value for this ratio is either a large value of nig or a large value for the 
time To; we prefer the former, and will express all other time and frequency scales in terms 
of tq. We make oc throughout the simulation, that is, we assume that the QCD 
scale (where stops varying) is reached after the simulation ends. We terminate each 
simulation when = 0.1, so the axion mass is coming on order the lattice spacing. At 

this time, we evaluate the axion-number content of the field and extract K. 

We begin with 2D simulations because they are numerically cheap, so we can achieve 
quite large values of m^ro. We look first at how the string and wall densities vary with 
time. Our results, scaling out the expected A'str oc , Awaii oc behavior, appear in 
Fig. 5. Much of the figure is as expected. The scaled string and wall densities rise with r 
below about r = 1.4ro, due to logarithmic scaling corrections we have already discussed. 
They are also larger for larger values of rrisTQ , also as a result of logarithmic corrections to 
scaling. The wall length starts to fall at about r = 1.2ro, and the strings begin to move 
faster around r = 1.5ro, leading to a brief peak in ^ as the strings’ 7 -factor rises and then 
a fall in the string length and extent starting around r = I.Dtq. Eventually, the extent of 
both strings and walls falls near zero. For larger values of mgTo, the strings are heavier 
by a factor of In(msro), and so they respond to the walls with more inertia. This delays 
slightly the growth in string velocity and extent, and slightly delays the fall of the wall and 
string extent. The 7 -factors we find at late times, while large, are actually underestimates, 
as the 7 -factor is especially sensitive to the lattice spacing, as we explain in Appendix A. 2 . 

But the figure shows something that might at first be unexpected. For each value of 
rUsTo, there is a point where the wall area abruptly crashes. At the same value, the string 
number briefly spikes, and then crashes as well. The larger the value of nigTo, the larger 
the To value where this occurs. In every case it occurs when rrialmg = 1/(43 ± 3). 

This collapse of the wall network is a numerical artifact which we will explain in 
Subsection 4.3. It occurs because the ratio m^/m^ is finite in these simulations, and when 
it becomes large enough, about = 1/39, the walls cease to be even metastable and 

become absolutely unstable. This leads to an abrupt unraveling or collapse of the network. 


- 13 - 



2D Wall Area, rrig case 




2D String speed, case 



Conformal time t/t^ 


Figure 5. String and wall density in 2D simulations, for a number of choices of msTQ. Top left: 
string number. Top right: wall length. Bottom left: string ^-factor. Bottom right: velocity, 
and 7 factor of strings. 


But this collapse has nothing to do with the physics that would occur in any circumstance 
where m^/m^ remains large. We should question the physical relevance of any simulation 
where this occurs before the string-wall network has largely broken up via the expected 
physics. We see that this constrains us to consider quite large values of m^ro. 

To see this a little better, we plot the axion number determined via Eq. (4.4) (full 
details in Appendix B) for all but the two largest nisTQ values in Fig. 6. The figure clearly 
shows that larger nigTo values, meaning denser string networks, lead to larger axion number 
at intermediate times. As the string network breaks up, the axion number falls. When the 
network abruptly unravels, the axion number drops sharply, and then goes flat after the 
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Figure 6 . Axion number from Eq. (4.4) versus time for simulations with msTo between 75 and 
900. 

network is gone. Except for the smallest nisTo, the final value is almost unchanged. Still, 
the figure leaves us distrustful of overinterpreting simulations where the string network 
crashes rather than breaking up via natural string-wall network dynamics. Roughly, the 
rUsTo = 900, 600, and 450 simulations look reasonable, but the abrupt drop is rather large 
for the smaller-msTo simulations. 

As a final way of looking at the network, we plot the energy content as a function 
of time in Fig. 7. The energy content naturally scales as /^/t^ on dimensional grounds 
and due to Hubble expansion, so we have factored this out. Also, at late times the rising 
mass causes the energy to rise, £ oc ma', but this behavior only sets in once the network 
starts responding to nia, at about r = 1.4ro. Therefore we have also scaled out a factor 
of (5 + maTo), with the 5 chosen so the scaling will turn on at about r = 1.4ro. To avoid 
crowding the plots we have not shown error bars. The largest m^ro curves have statistical 
errors around 2%; the smaller m^ro have smaller statistical errors. In the energy-fractions 
plot we have shown only the three largest m^ro values. The string energy is estimated 
as In(msr) with -|- r“^, so the IR cutoff is either the inverse string 

separation or the axion mass, whichever is larger. The domain wall energy is estimated 
from domain wall area, without an attempt to include the y-factor for the walls’ motion, 
thus the wall energy is probably an underestimate. The string energy is an underestimate 
after 7 gets large around r = I.Stq, since our method underestimates large string velocities 
when using such large m^a values. 

The figure shows that the energy starts out strongly tUsTo dependent, presumably 
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Scaled Energy, 2D simulations 


Scaled Potential Energy, 2D simulations 



Scaled Gradient Energy, 2D simulations 





Conformal time t/tq (mg(TQ)TQ=l) 


Figure 7. Energy density, scaled by t^/(/^(5 + niaT)), for 2D simulations. Top left: total energy. 
Top right: potential. Bottom left: gradient energy. Bottom right: energy fractions and estimated 
energy in strings and walls. 


because of the different density and string tension of the string network. By the end of 
the simulation this m^ro dependence has largely disappeared. Also, the early behavior is 
dominated by gradient energy, since strings at rest carry almost all their energy in gradients 
(moving strings have a kinetic energy fraction of n^/2). The late behavior is oscillating 
axions, with energy equipartitioned between kinetic and (potential+gradient). Quadratic 
fluctuations have 1/2 their energy as kinetic, whereas walls and strings have most energy 
in gradients and potential; so the extent that the (phase-averaged) kinetic energy is below 
1/2 is a reasonable estimate of how much energy is still in strings an walls. Once the energy 
is mostly in quadratic fluctuations (particles), there is still a shift between gradient and 
potential energy, due to the increase in rua- The oscillations at late time show that the 
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produced axions are not a phase-random collection, but have some phase coherence. 

We should worry about the impact of the abrupt collapse of the string-wall network if 
the strings and walls still carry significant energy fraction when it happens. The bottom 
right figure gives a nice criterion for knowing if this happens. In each case shown in the 
plot, less than 5% of energy in walls when they collapse. This is not true for smaller rngTo 
values, so the associated simulations cannot be completely trusted. 

4.3 Domain-wall instability 

The most striking feature of the string-wall network simulations we just presented is the 
sudden collapse in the amount of domain wall, accompanied by a brief spike, and then 
collapse, in the amount of string. The time value t/tq when this occurs changes as we 
vary rngTo. But it always occurs at the same value of m^/m^ ~ 1/43. Here we show that 
this is not an accident; it happens because the domain walls are only metastable, and at 
m^/m^ = 1/39 they become completely unstable. 

A domain wall is where the phase 6a changes from 0 to 27r by going around the 
“valley” of the winebottle potential. Since 9a is a single real parameter, this occurs on a 
codimension-1 surface. The thickness of the region where 6a is far from 0 or 27r is set by 
a competition between potential energy, which wants to make the region thin to keep the 
region with potential energy costs small, and gradient energy, which wants to make the 
region thick to minimize f The thickness of the region is ~ l/nia', when maT ^ 1, 

this is thin compared to the horizon scale and the walls can be approximated as planar. 

In the case <C m^, we can take \/2(p = fae^^°' to good approximation. Choosing 
the 0 a-variation along the 2 : axis, the energy to minimize is 

E^aii = j dxdy j dz {^{d^daf + x(l - cos^a)^ (4.6) 

= ^ = /f / dz (^^{d,6af + mlil - cos 6a)^ , 

with boundary conditions da -^z^-oc 0, 6a —>-z^oo 27r. Here A is the area of wall we 
consider, and a = E/A is the surface tension. Extremization gives 

6a = rua sin 6a = , (4.7) 

with V = 771^(1 — cos0a) = y/fa- This can be solved by multiplying both sides by dz6a, 

dz6adl6a = dz6ade,V{6a), (4.8) 

. {dz6a? d6a dv dv 

^ 2 dz d6a dz ’ 

and integrating f dz to obtain a Virial-type relation 

]j^{dz6af = V(6a) -V(0) = ma(l - cos 6a) , (4.9) 
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showing that potential and gradient terms each represent half the wall’s energy. The surface 
tension is® 


a 

~2 


d6a 


f2 
J a 


' —oo 
r2'K 


dzid.Oaf = / dz^^2V-2Vo 


p27T j - ^27r 

/ d9a\j2V{9a) - 2V{0) = rUa d9a-\/2 - 2 cos 9a = SrUa- (4.10) 

Jo Jo 


Numerically, we estimate the energy in domain walls as Siriafa times the area where 9a = vr, 
which we determine as described in Appendix A. 

Now suppose that m^/m^ is large but not enormous. Then (p{z) will not strictly lie 
on the circle /ae*®“, and Eq. (4.6) becomes (introducing ip = (f/fa) 


f2 
J a 

V{(pr,(pi) 


j dz Q {[dzifrf + [dz'fif) +V{(pr,(pi)'^ , 
-j) +”^a(l-^r). 


The equations of motion are 



dV 

diprp 


As z is varied, the field (p will follow a curve in the complex p plane. 


(4.11) 

(4.12) 


(4.13) 


p{z) = p{J{z)), 


(4.14) 


where i is the affine parameter describing the curve that (p follows; \dp/di\ = 1 and 
\dz(p\ = dl/dz. The complex equation of motion, Eq. (4.13), can be decomposed into the 
(r, i) component tangent to the curve and the component normal to the curve. The tangent 
EOM reads 

dll = d,V {dd? = 2{Vm))-%). (4.15) 

integrating into a Virial relation in the same way as before. The surface tension is again 

= I d£^2V{m)-^V{0). (4.16) 

The equation of motion in the field direction normal to the curve is what determines 
the curve ip{£); it reads 

(a _ 9V{p) d‘^p{i) _ dV{p)/d(pn , . 

dl^ dpn d£^ 2 V{p{£))- 2 Vq' ^ ’ 

This equation is equivalent to asking: what curve (p{l) will produce the lowest surface 
tension in Eq. (4.16)? By choosing a path with a small value of V, we keep the integrand 

^Some references use the value 9.32 rather than 8. This estimate originates from [58], who find it for 
the domain wall at T = 0, essentially by incorporating corrections to the strict (1 — cos 0a) form of the 
potential. Such corrections arise at zero temperature because instantons form a correlated liquid; but at 
high temperatures the dilute instanton gas approximation should be good and we expect such corrections 
to be tiny. 
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small; but by choosing a short path, we keep the integration range short. The LHS of 
Eq. (4.17), d‘^(p{l)lde, is the extrinsic curvature of the path. The larger its value, the 
more we can shorten the curve by rounding it off at this point. This must be balanced 
against how fast V will rise when rounding off the curve. The equation tells the exact 
balance between the gain of shortening the curve and the cost of increasing V along the 
curve. Now V—Vq in the denominator is mostly provided by m^, while the normal derivative 
in the numerator is almost purely provided by The larger we make the less 

there is to stop the domain wall from rounding off its path and exploring (pi + (pf < 1. 


Path of domain wall through ^-space 



¥’r/fa 


0.25 

0.20 



Ri- v>\/2//« 


Figure 8. Left; domain wall solution, as seen in (^-space, for several values of m^/rUg. Each curve 
is the path (p follows through field-space for a given mUml value. Neighboring dots are points that 
are separated by l/( 4 ms) in coordinate space. Right: the largest-m^/mg path as it appears on the 
(pr,Pi) potential, illustrating how the curve departs from the “valley” of lowest potential. 


We have solved explicitly for the domain wall’s shape and surface tension, using V 
from Eq. (4.12) with different values of We illustrate the results in Fig. 8. The 

left plot in the figure shows several (p{£) curves in the (pr, p>i plane; the curves shown have 
ml/ml values spaced in intervals of 0.0032. The dots are the values of (p{l{z)) at a series 
of z-coordinates, with neighboring dots on a curve separated by Az = l/(4ms). The wall 
becomes thinner as the dots become more spread out. The last curve is the last metastable 
domain wall; a tiny further increase in ml will cause the curve to pull over the top of the 
potential peak, and the domain wall spontaneously collapses. On the right in the figure, 
we have shown the path of this domain wall in terms of the potential V{(pr, p>i)^ indicating 
how the curve “pulls up” partway onto the bump in the potential rather than staying in 
the “valley” of the potential. 

Our study finds that the domain walls lose their metastability when ml/ml reaches 
about 1/39. For larger ml/ml values, there is no longer any domain wall solution. In 
simulations we observe the collapse of the wall area to start when the ratio is ~ 1/43. We 
believe that this is because fluctuations in the field, hitting the wall, can induce a collapse 
when ml has not quite reached the value where instability occurs. This will not happen 
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everywhere at once, but only locally where larger fluctuations impact the wall. At these 
spots the wall will “break”; a loop of string forms the boundary of this break, temporarily 
raising the total length of string in the simulation. The hole in the wall then rapidly grows 
and is joined by new breaks, leading to the collapse in the wall area; as the holes in the 
wall grow and percolate, the amount of string at first rises but then falls essentially to zero 
as the wall network disappears. This is a good description of both the timing and behavior 
of the wall collapse we observe. 

Physically, small values of mg are experimentally excluded. And it is most natural to 
expect A ~ 1 so mg ~ fa, which should be orders of magnitude larger than rua- So while 
this physics is the correct dynamics of a string-wall network with m^/m^ ~ 1/40, it does 
not describe the evolution of physical interest. Therefore we cannot rely on the results 
of any simulation in which a significant amount of energy still resides in the string-wall 
network when this collapse occurs. This criterion places a limit on what values of rUgTQ give 
reliable answers, pushing us towards fine lattice spacings and towards the largest values of 
niga that still give continuum behavior. 

4.4 3D simulations 

Our experience with 2D simulations tells us that we must consider the largest possible 
values of nigTo to avoid the unphysical collapse of the string network. This means we need 
to choose the largest m^a we can, subject to the constraint m^a < 1.5, found in Sec. A.3 to 
ensure continuum behavior. We should also choose the largest aro value we can, subject to 
the constraints that t/tq must get large enough for the network evolution to complete, with 
L > 2r to ensure no finite-volume errors. Unfortunately, our limited numerical resources 
limit us to boxes of 1600^ or smaller, constraining us to consider aro = 300 or nigTo = 450 
but not larger. We have not considered m^ro smaller than 225, since the 2D simulations 
showed that the wall network breaks up too early in such simulations to learn anything of 
value. So we have less dynamic range than in the 2D simulations. 

The 3D simulations show surprisingly similar behavior to those in 2D. Rather than 
repeating all plots we presented in 2D, we plot just the string length and speed, the wall 
area, and the energy ratios in Fig. 9. The most signihcant difference from 2D is that, 
while the strings and walls start to decay at about the same time, the network decays more 
quickly, with very little string left by r = 2.57ro, when the network collapse occurs for 
nigTo = 450. In particular, the collapse of the walls, in the plot of wall area and of energy 
ratios, is almost invisible for this largest nigTo value. Therefore, while this final simulation 
does not provide the right string tension due to the missing core tension, it at least presents 
a case where the network breaks up via a physical mechanism rather than the loss of wall 
stability. 

For further comparison between the 2D and 3D cases, we have plotted the energy and 
the energy fractions for each dimensionality at a common value of nigTo = 450, in Fig. 10. 
The energy density starts out somewhat higher in the 3D case, corresponding to the larger 
string density obtained in 3D; but the later stages of the evolution are strikingly similar, 
except that the 3D string/wall network decays before the unphysical collapse, while the 
2D network still carries some energy when the collapse occurs. 
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Energy fractions, 3D (am,= 1.5 lgj^^^j.=2.25/mQ) 





Figure 9. String length, wall area, string velocity, and energy ratios as a function of time in 3 
dimensions for 3 values of TOsTo- 

We also plot the final axion-production efficiency K as a function of nigTo, for both 
2D and 3D, in Fig. 11. Besides the very smallest values of m^ro, it is a very weak function. 
Surprisingly, K is about half of the angle-averaged misalignment-mechanism estimate, and 
it appears to be a weakly declining function of m^ro. Note however that it is very dangerous 
to extrapolate based on this result, since the physical value of rrisTo 10^'^ is about 27 
orders of magnitude larger. 

5 Discussion and Conclusions 

The axion abundance is determined by dynamics occurring around the conformal time tq 
where Tm{T) = 1; we have seen that the axion number is essentially fixed by r = Stq. 
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Scaled energy components, 2D vs 3D 



Conformal time t/tq 


Energy fractions. 3D vs 2D, m,TQ=450 



Figure 10. Comparing 2D and 3D simulations at the same uIsTo value: Scaled energy (left) and 
energy fractions (right). 
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Figure 11. Axion production efficiency K as a function of itIsTo, in 2D and in 3D 

Entering this time period, there is a network of axionic cosmic strings, which evolves and 
breaks up due to the appearance of the axionic mass. 

We have shown that in both 2D and in 3D, the initial density of the cosmic string 
network is a strong, roughly linear, function of the log of the horizon-to-core ratio In(msro). 
In simulations we can make ln(m 5 To) as large as ln(450) = 6 in 3D and ln(1800) = 7.5 in 
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2D; but in nature we expect it to be ln( lO^O) rsj 69. Therefore we can only study axion 
production (so far) for cases with a much more dilute starting string network than we 
expect in the cosmological axion context. 

We also showed that, for too small a value of m^ro, the network dynamics suffer an 
additional numerical artifact; the axionic walls and strings become absolutely unstable and 
the network suddenly collapses. This occurs whenever the ratio of angular-to-radial masses 
exceeds m^/m^ = 1/39. This effect compels us to work at large (numerically expensive) 
values of rngTo; but we can achieve values large enough that the instability only occurs 
after almost all strings and walls have disappeared. 

We do not believe that this problem influenced the results of [17] because their m^/m^ 
ratio stopped increasing partway through their simulations and then remained fixed, pre¬ 
venting it from growing too large. On the other hand, for the axion to play a part in dark 
matter we expect tq to occur at temperature Tq ~ 1.5 GeV, so this flattening-off behavior 
is unphysical, raising some questions about their results. 

Perhaps the most significant, and to us surprising, result of our simulations is that, 
even over a range of k = In(msro) where the starting string density varies by a factor of 2, 
the final axion number is unchanged at the 20% level, and in fact trends slightly down with 
increasing string network density. So at least within the range of string network densities 
we have been able to study, the string density appears to play little role in establishing the 
final axion number density. 

Suppose this is the case, so the produced axion number is well described by Eq. (4.3) 
with K = 8. What do we learn about the axion decay-constant fa and mass m^? Combin¬ 
ing Eq. (4.3), the Friedmann equation for a hot gas. 


Stt 

3mpj 30 


(5.1) 


an estimate of the hot topological susceptibility from Wantz and Shellard [46], 


X{T » To) ~ a^s^\K/TT , 


(5.2) 


with A = 400 MeV, n ~ 6.68, and Ows = 1-68 x 10 and the relevant cosmological 
information from Planck [38] 

— ~ 8.59 X 10"^^ , 
s 


Pdm 

s 


DM ^ 

s 


^'^^^^ (938 MeV)(8.59 x 10"^M ~ 0.39 eV, 
0 . 0221 ^ ^ 


(5.3) 


along with an estimate g^, = 64 based on a plasma of {e/iT, i'i 23 ,'y,g,udsc) with the QCD 
degrees of freedom contributing 80% of the free-particle value to account for strong inter¬ 
action corrections at this temperature [59], we find 


fa = 3.3 X 10^^ GeV x 


\ 8 ) V64y Vl-68xl0-7 


0.079 


(5.4) 


This value gives rua = (76 MeV)^//a = 18/reV. The transition temperature is Tq ~ 1.5 
GeV. In our simulations, most of the nontrivial 3D network evolution took place between 
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r = 1.4to and r = 2.6to, corresponding to T = 1100-580 MeV. This is the range where we 
need to understand the topological susceptibility better - though fa is only sensitive to a 
rescaling of the topological susceptibility through the 0.079 power, so even a factor of 10 
error in the estimate of Eq. (5.2) (as suggested by recent work [48]) makes a modest 20% 
shift in fa and nia- A significant change in our estimate for K would have a larger effect 
(though a very small effect on the relevant Tq). 

Over the range we studied, there was very little change to K in going from 2D to 
3D and in varying the string network density by about a factor of 2. We have made no 
attempt to separate which axions arise from misalignment, which from strings, and which 
from walls; indeed it is not clear to us that doing so is either well defined or terribly useful. 
But it appears that there is not a large component strictly proportional to the density of 
strings. Nevertheless, we find it rather brave to extrapolate that the same value K = 8 
should apply when the string network is 5 to 10 times denser than in the simulations we 
considered. 

We believe that it is well motivated to look for a way of simulating axion production 
from global networks with much larger core tension. It is clear that the enormous factor we 
need cannot be achieved simply by shrinking the lattice spacing. Rather, a new strategy 
is needed. We see the need for a method to excise and treat explicitly the physics of the 
string core, by adding degrees of freedom to describe its tension and inertia. We are close 
to presenting such a technique for 2 dimensional networks, taking advantage of the dual 
electromagnetic description. Details will appear elsewhere. 
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A Nnmerical implementation 
A.l Action approach 

The starting point of our implementation is the action, Eq. (3.2), with —)• in the first 
potential term. We will directly discretize the action; the extremization with respect to 
each field value then gives us an update rule which will be a leapfrog update. It is also 
straightforward in this approach to make the temporal spacing vary, for instance having 
finer time step at very early times when the behavior (which gives Hubble drag) is 
rapidly changing from timestep to timestep. Specifically, if we discretize on a cubic lattice 
with spacing a and a set of conformal times with spacing Or,n = — Tn-i, then the 
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action is 


C< _ 3 / TnTn—l 

^ ^ ® 0‘T,n 


((Prix,Tn) - ipr{x,Tn-l))‘^ + (r z) 


(A.l) 


+EE 
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+EE 


3 2 ^T,n+1 


a r. 
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,i=l,2,3 L 


(ipr(x + ai, Tn) - <frix, Tn)f + (r z) 


3 2 ^r,n+l / r 2 


a T, 


[ipj{x,Tn) + ipUx,Tn) - l] + m^(T)(l - ifr) 


Here the time derivative term stretches between times Tn-i and t„, so we replace in its 
coefficient with the product of the starting and ending time. The coefficient ar,n on this 
term accounts for the f dr running from r^-i to Tn- For the space terms, which occur at 
time Tn, we replace with r^, and we treat the amount of f dT contributing to the term 
to be half the interval before the term appears plus half the interval after the term appears, 
hence the factor {a-r^n + ar,n+i)/2. We have also implemented a next-neighbor improved 
gradient term, but most of our results are based on the above gradient term. 

For uniform temporal spacing, variation with respect to y?, gives 


‘-Pi{x,Tn+l) - ipi{x,Tn) 


IH-ILA Tn) - (Pi{x, Tn-l)) 

'XnXn+l 


+ - 


XnXn+1 



Piix + ai,Tn) - 2<pi{x,Tn) + Pi{x - ai,Tn) 


a2 


"is 2 

+ “ ‘Pr 



(A.2) 


One often writes this in terms of the conjugate momentum ipi{x,Tn+i) — ipi{x,Tn) = 
'Ki{x,Tn)- The factor Tn-i/Tn+i, which arises from the overall factor in the action, 
accounts for Hubble drag. The middle term is the gradient term, the final term is the 
radial potential term, and the (pr equation is the same but with the addition of a linear 
term. 

In practice we use a very fine temporal spacing at small r, a,- = a/40 for r < a and 
a-r < r/20 thereafter, to avoid errors when the Hubble drag is large. This is probably 
unnecessary. At larger times we go over to a fixed value for ar/a. 

As initial values, we first set the field to be of unit magnitude and independent random 
phase {ipr{x) = cos^a, and Pi{x) = sin^a,, with each Ox chosen uniformly from [0,27r)) 
and then apply smearing to a coherence length fsmear- We will check for dependence on 
parameters such as ztZsO, 4mear, and at/a momentarily. 

A.2 Strings and string velocities 

We want to know the density of the string network, which requires identifying where the 
strings are. To do this we define the plaquettes that are pierced by a string, and we count 
these plaquettes, applying a statistical correction which we now explain. 
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Each plaquette has four corners, and we say that a string goes through a plaquette if 
the tetragon in y? space, whose corners are the if values at the corners of the plaquette, 
encloses the point (^ = 0, see Figure 12. The direction or sense of the string is determined 
by whether the origin is enclosed in a clockwise or a counterclockwise sense. 



Figure 12. Four points around a plaquette map to four points on the complex ip plane. If the 
resulting tetragon encloses the origin, we identify the plaquette as being pierced by a string. 

Our algorithm for determining if a string pierces a plaquette is as follows. For the 
tetragon in Fig. 12 to contain a string, the real axis must be crossed twice with the same 
handedness (clockwise or counterclockwise); that is, the number of strings piercing the 
plaquette is half the signed sum of real axis crossings, where the sign is determined by 
whether the crossing is clockwise or counterclockwise. The real axis is crossed between 
points xi and X 2 with values fi = ip{xi) and ip 2 = f{x 2 ) if Ini (<^ 1 ) Im (<^ 2 ) < 0. The 
axis crossing is clockwise if \m{ipiip 2 ) > 0 and is counterclockwise if Im((/9i(^2) < 0- The 
winding number is the sum of (~|) each clockwise (counterclockwise) axis crossing, 
as we consider each pair of corners going clockwise around the plaquette. 

We also identify which links have a domain wall go through them by finding links 
where Im ((^Ji) Im ( 992 ) < 0 (the real axis is crossed) with lm{ipiip'^lm.{ipi — (^ 2 ) < 0 (it 
is the negative half-axis which is crossed). This occurs an odd number of times for each 
plaquette containing a string, and an even number of times for any other plaquette. Each 
algorithm (string and wall) requires only multiplication and comparison; no divisions or 
trigonometric functions are needed, providing good numerical efficiency. We can sample 
the full lattice for strings every lattice unit of time (Ar = a) while taking less numerical 
effort than the field updates. We believe that our plaquette identihcation is equivalent 
to that of Hiramatsu et aJ [16], but by avoiding explicit reference to angles we avoid the 
trigonometric evaluations needed there. 

We make no attempt to connect plaquettes pierced by strings to identify the strings 
and their lengths. Instead we rely on counting plaquettes with string and links with wall. 
In 2 dimensions this works for strings, but in 3d (and in both 2d and 3d for walls) there 
are normalization issues. If the unit normal along a string is then a length 

L of string will pierce L\lx\ yz-plaquettes, L\ly\ xz-plaquettes, and L\lz\ xy-plaquettes; 
for different directions the total number of plaquettes is between L and Ly/S. We are 
only trying to determine the string density statistically, and all directions are statistically 
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equally likely, so we can find the right average string density statistically by counting 
string-pierced plaquettes and dividing by the angle-averaged number of plaquettes per unit 
length of string, which is 

count-factor = [ — [ — -= 3 [ dcos{9) x cos{9) = - , (A.3) 

Jo 27r Jq 2 Jq 2 

where ^x,y,z are the x, y, z-components of the unit vector in the 4>, 9 direction. The same 
overcounting factor applies for domain walls in 3D, since the number of links piercing the 
wall depends on the wall-normal in the same way as plaquettes depend on the string-normal. 
In 2D the domain walls are overcounted by -|- |Dy|) = d/vr. 

To determine the average string velocity v, and 7 = (1 — we use the time 

derivative of the field near the string core. Consider a string core stretched along the z axis 
and moving in the x direction. For the string at rest, the minimum of the string’s energy. 


Eq. (2.1), occurs for 

ip{x, y, z) = fa^-^f{r) , (A.4) 

with /(r) solving the equation of motion 

^^-^/(l —/^) = 0 , /(0 ) = 0, /(r — 00 ) —)• 1. (A.5) 

The equation of motion near 0 enforces that 

r/\ \S C“l“ 16c^ / / A \ 

J [r) = cmsT -—[msr) -\ -—(m^r) -h ... (A.6) 

with c a constant; solving Eq. (A.5) via overshoot-undershoot determines c = 0.41238. 
Therefore, near the string core, the field is 

y) = facrusix + iy)e~''^° |^1 - .. .^ . (A.7) 


The moving string solution at t = 0 is found by replacing x with 7 a:; and the time derivative 
for the moving string is dt^p = —vdxP- Therefore 


iflin o / N ^ f 3 f 37 ^x 2 -h 2i'yxy + y‘^\ 

- e y\v) = jacmsjv - facm'l'yv I - — - I 

At lowest order in mJr'^ we hnd 

o /?2 2 2 2 2 2 2 

dtp dtp = fam^c j V 7 u = 


-F . . . . 


dtp*dtp 


(A. 8 ) 


(A.9) 


The next-order term in Eq. (A. 8 ) can be used to engineer a correction for fields close to 
but not at the center of the string: 


2 2 _ dtp*dtp f p*p \ {p*dtp + pdtp*f 

~ V 8c^fa) 


(A.IO) 
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which we found by establishing the leading-order small-distance behavior of ^p*(p and 
\ip*dt^\ and finding a combination that would cancel the subleading contributions in 
dt(p*dt^. 

We apply this approach by sampling over points that are near the string’s center, so 
that the second-order, ~ terms are small, and higher (uncomputed) corrections should 
be negligible. We choose for our sample of points near the string core the set of points on 
corners of plaquettes that are pierced by a string; points on more than one string-pierced 
plaquette are counted once per plaquette. For each plaquette we find u^/(l — explicitly 
by averaging Eq. (A. 10) over the four plaquette corners; we then compute the local value 
of (v), (u^), and ( 7 ). We average (l,u,u^, 7 ) over all pierced plaquettes, weighting with 
weight 7 in order to follow the literature convention that the string network should be 
weighted by f 'ydl (energy content), not f dl (length). The numerical overhead is small 
because the computation need only be performed on sites that are identified as corners of 
a string-pierced plaquette. 

The method assumes that points with distance ~ a away from the string core are still 
“relatively near” the string’s core, which is an expansion in (cmsa)^/16 < 1. Note that, 
in 3D, the average distance-squared of a point on a pierced-plaquette to the nearest point 
on the piercing string is (x^) = a^/ 2 , so the sampled points are surprisingly close to the 
strings, and the higher-order corrections are expected to be small. 

Another weakness of the approach is that it is based on the structure of a straight 
string; it will commit errors for bent, accelerating strings. However, in this case the string’s 
length and velocity are anyways not uniquely defined; it is not clear to us that the method 
works any worse than other approaches. Similarly, it could be confused by any radial 
fluctuations (breathing modes) in /(r); but these modes are heavy and we do not expect 
them to occur with large amplitude except perhaps at first due to initial conditions. Also 
note that, since the algorithm determines in terms of a manifestly positive expression, 
it never concludes that a string is moving faster than light speed, something which could 
in principle occur for some string-velocity approaches. 

We mention one final challenge for our method. We have assumed that the string’s 
shape correctly reflects the Lorentz contraction of the string’s structure. On a sufficiently 
coarse lattice, Lorentz invariance is not respected in the string core and the Lorentz con¬ 
traction will not be as strong as it should be in the continuum. This means that the method 
will under-estimate the 7 factor of very fast strings on coarse lattices. If the main goal is 
a high-precision determination of string velocities and 7 -factors, one should use a smaller 
value of mstt than we have used in most of our studies. 

To test this last point, we have repeated the simulation of a nisTQ = 450 2D string 
network using m^a = 0.75 and mga = 1.0 lattices, to compare with our results in the 
main text using mga = 1.5. Most properties (energy, string length, wall extent, final axion 
number) agree very well; the final axion number is the same within 1 % statistical errors. 
But once the strings start to move very fast late in the simulation, the finer lattice observes 
a significantly larger 7 -factor for the strings, as shown in Fig. 13. 

We think a similar velocity-finding approach could be applied to field-theory simula¬ 
tions of local strings (Abelian Higgs model simulations), but since the string’s structure 
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Figure 13. String velocity for 2D, TOsTo = 450 simulations with three values of lattice spacing: 
TTisO = 0.75, 1.0, 1.5. Our string velocity method finds a higher y-factor on the hner lattice. 
(The crash in the y-factor at r = 2.57to coincides with the collapse of the network by direct wall 
instabilities and should not be taken seriously.) 

then depends on two parameters (the Higgs mass and the gauge boson mass), the applica¬ 
tion is more complicated, with the constant c above replaced by some dependent 

value; finding the NLO corrections would also be significantly more complicated. 

A.3 Numerical tests 

We want to find the largest values of mga and at/a that are compatible with a continuum 
interpretation, and we want to check for sensitivity to initial conditions such as ^smear- We 
will test these parameters using ttIq = 0 or string-only simulations, because there are then 
fewer scales to consider. 

We start with m^a. If we choose this too large, the string core is < 1 lattice site 
across, which is too small for the lattice to resolve properly. In this case, the UV edge 
of the integral in Eq. (2.2) becomes sensitive to the exact location of the string relative 
to the lattice, leading to a string energy varying periodically with period a, rather than 
being translation independent. This makes it possible for a string to “stick” in the most 
energetically favorable lattice location, which will interfere with string evolution and impede 
the annihilation of the network. To test for this problem, we made a series of evolutions 
with different values of mga, but identical other properties as measured in terms of mg- 
These can be interpreted as varying the lattice spacing. 

We show how the string length (or, in 2D, the number of strings) vary with time for 
several values of m^a, in both 2 and 3 dimensions, in Fig. 14. In both cases, the results are 
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Figure 14. Dependence of string density on the lattice spacing (in units of the saxion mass, amg). 
Left: string length in 3 dimensions. Right: number of strings in 2 dimensions. Lines represent Icr 
statistical range based on averaging several samples. 


consistent within 2% for all mga < 1.5, but deviate for larger m^a, especially at later times. 
We did a similar study with the next-neighbor improved action (not shown), which showed 
continuum behavior slightly sooner, at msO, = 1.8. The reduced number of lattice points 
this allows, (1.8/1.5)'^"''^ in d space dimensions, does not quite make up for the factor of 
2 in numerical cost for that algorithm, so we have generally stuck with nearest-neighbor 
interactions in the remainder of our work. 

Note that the result in Fig. 14 clearly shows that the string length rises as In(mst), 
rather than approaching a flat value, both in 2D and in 3D. Therefore we already see the 
logarithmic corrections to scaling in this figure. If we re-plot the figure in terms of r/a, 
rather than rUsT, the lines do not fall on top of each other; the relevant physical length 
scale for comparison is mg, which sets the string core size, not the lattice spacing. 

Next we check for independence on initial conditions, by varying the length scale ^smear 
over which the initial conditions are smeared. The result, shown in Fig. 15, indicates that 
different initial conditions rather quickly converge to the same string network density, which 
again scales logarithmically with the system age as measured in units. This means the 
choice for this smearing length is not important. We typically choose 2.1/mg in this work. 

We should also check what temporal spacing is sufficient to approximate continuous 
time. There is a Courant condition beyond which the update algorithm becomes unstable: 
Of/a = 2/y4d + mJ with d = 2,3 the number of space dimensions. We choose a/at integer, 
so this limits the available values to a/at = 2,3,.... The axion is very light, and at late 
(interesting) times we expect that the physics is primarily in long wave lengths except near 
string cores; so one might expect less severe dependence on this ratio than for some lattice 
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Convergence to Scaling Behavior 



Figure 15. Dependence of string density, in 2D, on the initial correlation length of the random 
initial conditions. 




Figure 16. Temporal spacing dependence of string density (left) and energy (right). In each case 
we compare with an extrapolation to zero spacing, based on Ot/a = 1/6 and 1/8 and assuming 
(ot/a)^ scaling of errors, which is the expected scaling form. 


problems. This turns out to be the case; as we show^'^ in Fig. 16, a temporal spacing of 

suppressed statistical fluctuations by using identical initial conditions for each at value we con¬ 
sidered. The noise in the lines on the right in the figure arise because we write out too few digits before 
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atja = 1 /2 actually only leads to percent differences in string length and energy, compared 
to small temporal spacing. Nevertheless, out of paranoia we have generally used at/a = 1/6 
or 1/8 in other measurements, which should keep at errors below the per-mille level. 



0 aOO 400 600 800 

Conformal time rm^ 

Figure 17. Volume dependence of string length, in 2 dimensions, indicating Icr error ranges; based 
on 800, 3000, 12800, 12800, 12800 samples. 

Finally, although there are strong theoretical arguments that the volume cannot affect 
the network’s statistical properties so long as r < L/2 (in conformal coordinates), we 
check this nevertheless in Fig. 17. The hgure shows that r must grow to nearly twice L/2 
before significant finite-volume corrections occur; but they then rapidly become severe. All 
results in the rest of the paper consistently use volumes large enough that evolutions stop 
at r < L/2. 

We have generally re-checked the 2D-only results of this section with 3D simulations, 
which give consistent results but with poorer statistics because the added numerical costs 
of treating the 3D system make it harder to achieve high statistics. 


B Counting axions 


At the end of a simulation we have from which we want to extract an axion number. 

To do so we first convert to the axion field amplitude 9a, and then extract the axion number 
stored in the field. The axion field amplitude 9a is determined as 

Im (p*dr^ 


9a = argy?, dr9a = 




(B.l) 


processing the data; they drop in size when the energy becomes one digit shorter so one more digit is written 
out. 


- 32 - 










Since the linear term, Eq. (3.3), shifts the absolute minimum of the potential slightly away 
from Lfr = fa, we make an addative shift to ipr so that its minimum is at fa before applying 
these rules. 

We then assume that this axion held obeys the quadratic Hamiltonian (writing 6 = drO) 
H{ea) = fa I (^6a{ml - V^)9a + 9^) . (B.2) 

The Hamiltonian is diagonal in Fourier space {V is the volume of space) 

^ = j + ^a(^)) • (B.3) 

The particle number associated with a mode with oscillation frequency oj = + rria is 

the energy over the frequency, so the particle number density is 


IT'Slx 


d^k 1 
(27r)3 2 


^fWT^9l{k) + 


1 

^k'^+ml 



(B.4) 


Evaluating this in Fourier space is straightforward and is efficient when the EFT is 
available. If the box size is not a power of two or if the data is divided over processors in an 
inconvenient way, we can use Laplace methods instead. It is easy to write an algorithm to 
evolve 9,9 in dissipative “time” f, with boundary conditions that 9{f = 0) is the original 
value of 9a: 

df9a{x,f) = - ml)9a{x,f) , 9a{x,f = 0) = 9a{x) , (B.5) 

and the same evolution for 9a (in terms of r evolution, 9a is considered an independent 
held). The reason to do so is that we can hnd 9a{x,f) by position-space methods, but we 
know that the evolution in Fourier space will be 


9a{k, f) = 9a{k) exp(-(A:^ + ml)f), 


(B.6) 


so large-A: modes are more quickly suppressed than small-A: modes. Note that 


[l n ^ ^ 1 

V TT Jo fk‘^ + ’ 


(B.7) 


so therefore 



(B.8) 


By evolving in r and numerically implementing the hrst integral, which can all be done in 
coordinate space, we get the desired second integral, which correctly counts axion number. 
We have compared this method to the FFT method where both are applicable and have 
conhrmed that they give the same answer up to controllable numerical (f-spacing and 
large-r extrapolation) issues, which are easily held below 1% with smaller numerical cost 
than the time evolution needed to find the final 9a configuration. However, this method is 
not practical for making frequent determinations of nax over the course of a simulation. 
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